!*****************************************************************************************
! This program contains the following functions for the estimation of the baseline model:
!  1. FUNCTION phi_update: compute the one-period-ahead belief
!  2. FUNCTION phi_bayes: update workers' prior
!  3. FUNCTION problevel: compute level probabilities (used in used in E_likegaussj)
!  4. FUNCTION seprate: exogenous exogenous probabilities
!*****************************************************************************************
MODULE mod_update


USE mod_types
USE mod_param_initial
USE mod_parameters
USE mod_classerror


IMPLICIT NONE


CONTAINS



!*************************************************************************************************
!(1) Function to compute the one-period-ahead belief (used in value iteration) with outcome = 1,2
!*************************************************************************************************
FUNCTION phi_update(phi,task,outcome)

REAL(dp), INTENT(IN):: phi
INTEGER, INTENT(IN):: task, outcome
REAL(dp):: phi_update



!------------
!a. Level 0
!------------
!firm's outside option (no updating)
IF (task==0) THEN 
phi_update = phi



!------------
!b. Level 1
!------------
!outcome 2 (high) 
ELSE IF (task==1 .and. outcome==2) THEN
phi_update = (1.0_dp - lambda_11)*alpha_1*phi/(alpha_1*phi+beta_1*(1.0_dp - phi)) + lambda_11*omega_21

!outcome 1 (low)
ELSE IF (task==1 .and. outcome==1) THEN
phi_update = (1.0_dp - lambda_11)*((1.0_dp - alpha_1)*phi)/((1.0_dp - alpha_1)*phi+(1.0_dp - beta_1)*(1.0_dp - phi)) &
               + lambda_11*omega_21



!-------------
!c. Level 2
!-------------
!outcome 2 (high)
ELSE IF (task==2 .and. outcome==2) THEN
phi_update = (1.0_dp - lambda_12)*alpha_2*phi/(alpha_2*phi+beta_2*(1.0_dp - phi))  + lambda_11*omega_22

!outcome 1 (low)
ELSE IF (task==2 .and. outcome==1) THEN
phi_update = (1.0_dp - lambda_12)*((1.0_dp - alpha_2)*phi)/((1.0_dp - alpha_2)*phi+(1.0_dp - beta_2)*(1.0_dp - phi))  &
              + lambda_12* omega_22



!-------------
!d. Level 3
!-------------
!outcome 2 (high)
ELSE IF (task>=3 .and. outcome==2) THEN
phi_update = (1.0_dp - lambda_13)*alpha_3*phi/(alpha_3*phi+beta_3*(1.0_dp - phi)) + lambda_11*omega_23

!outcome 1 (low)
ELSE IF (task>=3 .and. outcome==1) THEN
phi_update = (1.0_dp - lambda_13)*((1.0_dp - alpha_3)*phi)/((1.0_dp - alpha_3)*phi+(1.0_dp - beta_3)*(1.0_dp - phi)) &
                + lambda_13*omega_23




ELSE 
PRINT*, "Error, no new value for the one-step-ahead posterior"

END IF 

!phi_update = phi

END FUNCTION phi_update


!***********************************************************
!(2) Functions to explore the effect of drifts in beliefs
!***********************************************************
FUNCTION hk1(level_hk)

INTEGER, INTENT(IN):: level_hk

REAL(dp):: hk1


IF (level_hk==1) THEN

hk1 = lambda_11


ELSE IF (level_hk==2) THEN
 

hk1 = lambda_12


ELSE IF (level_hk>=3) THEN
 

hk1 = lambda_13

ELSE 

!PRINT*, 'error in calculating the belief drift'

END IF



END FUNCTION



FUNCTION hk2(level_hk, outcome, tenuretime)

INTEGER, INTENT(IN):: level_hk, outcome, tenuretime

REAL(dp):: hk2


IF (level_hk==1) THEN
 

!omega_11 = 0.0_dp 
!omega_12 = 0.0_dp
!omega_13 = 0.0_dp

!hk2 = lambda_11*omega_21*REAL(1-outcome,dp)
hk2 = lambda_11*(omega_21+omega_11*tenuretime)


ELSE IF (level_hk==2) THEN
 

hk2 = lambda_12*(omega_22 + omega_12*tenuretime)


ELSE IF (level_hk>=3) THEN
 

hk2 = lambda_13*(omega_23+  omega_13*tenuretime)

ELSE 

!PRINT*, 'error in calculating the belief drift'

END IF



END FUNCTION




!****************************************************************
!(3) Function to update beliefs used in the likelihood function
!****************************************************************
FUNCTION phi_bayes(phi,lev_hist,task1,outcome1,task2,outcome2,task3,outcome3,task4,outcome4,&
         task5,outcome5,task6,outcome6,task7,outcome7,task8,outcome8,task9,outcome9, priortype,ttenure)


REAL(dp), INTENT(IN):: phi
INTEGER, INTENT(IN):: lev_hist, task1, outcome1, task2, outcome2, task3, outcome3, task4, &
            outcome4, task5, outcome5, task6, outcome6, task7, outcome7, task8, outcome8, &
			task9, outcome9, priortype, ttenure
REAL(dp):: phi_bayes
REAL(dp), PARAMETER:: tinybayes = 1.0d-10

!-------------------
!a. Period-2 update
!-------------------
IF (lev_hist==1) THEN
phi_bayes = (1.0_dp - hk1(task1))*(  p_success(2,task1,outcome1, priortype)*phi  )&
           /(  p_success(2,task1,outcome1, priortype)*phi &
		      +p_success(1,task1,outcome1, priortype)*(1.0_dp-phi)  )  + hk2(task1,outcome1,ttenure)


IF (phi_bayes < 0.0_dp) THEN

 phi_bayes = 1.0d-10

ELSE IF (phi_bayes > 1.0_dp) THEN

 phi_bayes = 1.0_dp - 1.0d-10

!PRINT*, "The posterior in t=2 is not interior"

END IF    


!-------------------
!b. Period-3 update
!-------------------
ELSE IF (lev_hist==2) THEN
phi_bayes =  (1.0_dp - hk1(task2) )*(  (  p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype) )*phi  )&
           /(  ( p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype) )*phi& 
			+p_success(1,task1,outcome1, priortype)*p_success(1,task2,outcome2, priortype)*(1.0_dp-phi)   ) &
			+ hk2(task2,outcome2,ttenure)



IF (phi_bayes < 0.0_dp) THEN

 phi_bayes = 1.0d-10

ELSE IF (phi_bayes > 1.0_dp) THEN

 phi_bayes = 1.0_dp - 1.0d-10

!PRINT*, "The posterior in t=3 is not interior"

END IF     


!-------------------
!c. Period-4 update
!-------------------
ELSE IF (lev_hist==3) THEN
phi_bayes =  (1.0_dp - hk1(task3) )*(  ( p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)&
				  *p_success(2,task3,outcome3, priortype))*phi )&
           /(  (  p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)&
				   *p_success(2,task3,outcome3, priortype))*phi&
			+p_success(1,task1,outcome1, priortype)*p_success(1,task2,outcome2, priortype)*p_success(1,task3,outcome3, priortype)&
			*(1.0_dp-phi)  )  + hk2(task3,outcome3,ttenure)


IF (phi_bayes < 0.0_dp) THEN

 phi_bayes = 1.0d-10

ELSE IF (phi_bayes > 1.0_dp) THEN

 phi_bayes = 1.0_dp - 1.0d-10

!PRINT*, "The posterior in t=4 is not interior"

END IF     


!-------------------
!d. Period-5 update
!-------------------
ELSE IF (lev_hist==4) THEN
phi_bayes =  (1.0_dp - hk1(task4))*(  (  p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)*p_success(2,task3,outcome3, priortype)&
            *p_success(2,task4,outcome4, priortype))*phi  )&
           /(  ( p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)&
				   *p_success(2,task3,outcome3, priortype)*p_success(2,task4,outcome4, priortype) )*phi&
			+p_success(1,task1,outcome1, priortype)*p_success(1,task2,outcome2, priortype)*p_success(1,task3,outcome3, priortype)&
			*p_success(1,task4,outcome4, priortype)*(1.0_dp-phi)  ) + hk2(task4,outcome4,ttenure)


IF (phi_bayes < 0.0_dp) THEN

 phi_bayes = 1.0d-10

ELSE IF (phi_bayes > 1.0_dp) THEN

 phi_bayes = 1.0_dp - 1.0d-10

!PRINT*, "The posterior in t=5 is not interior"

END IF     
  


!-------------------
!e. Period-6 update
!-------------------
ELSE IF (lev_hist==5) THEN
phi_bayes =  (1.0_dp - hk1(task5))*(  (  p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)*p_success(2,task3,outcome3, priortype)&
            *p_success(2,task4,outcome4, priortype)*p_success(2,task5,outcome5, priortype))*phi)&
           /(  (   p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)*p_success(2,task3,outcome3, priortype)&
            *p_success(2,task4,outcome4, priortype)*p_success(2,task5,outcome5, priortype))*phi&
			+p_success(1,task1,outcome1, priortype)*p_success(1,task2,outcome2, priortype)*p_success(1,task3,outcome3, priortype)&
			*p_success(1,task4,outcome4, priortype)*p_success(1,task5,outcome5, priortype)*(1.0_dp-phi) ) &
			+ hk2(task5,outcome5,ttenure)


IF (phi_bayes < 0.0_dp) THEN

 phi_bayes = 1.0d-10

ELSE IF (phi_bayes > 1.0_dp) THEN

 phi_bayes = 1.0_dp - 1.0d-10

!PRINT*, "The posterior in t=6 is not interior"

END IF     
   

!-------------------
!f. Period-7 update
!-------------------
ELSE IF (lev_hist==6) THEN
phi_bayes = (1.0_dp - hk1(task6) )*( (  p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)*p_success(2,task3,outcome3, priortype)&
            *p_success(2,task4,outcome4, priortype)*p_success(2,task5,outcome5, priortype)&
			*p_success(2,task6,outcome6, priortype))*phi)&
           /( ( p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)*p_success(2,task3,outcome3, priortype)&
            *p_success(2,task4,outcome4, priortype)*p_success(2,task5,outcome5, priortype)&
			*p_success(2,task6,outcome6, priortype))*phi&
			+p_success(1,task1,outcome1, priortype)*p_success(1,task2,outcome2, priortype)*p_success(1,task3,outcome3, priortype)&
			*p_success(1,task4,outcome4, priortype)*p_success(1,task5,outcome5, priortype)*p_success(1,task6,outcome6, priortype)&
			*(1.0_dp-phi)  ) + hk2(task6,outcome6,ttenure)


IF (phi_bayes < 0.0_dp) THEN

 phi_bayes = 1.0d-10

ELSE IF (phi_bayes > 1.0_dp) THEN

 phi_bayes = 1.0_dp - 1.0d-10

!PRINT*, "The posterior in t=7 is not interior"

END IF     
   

!-------------------
!g. Period-8 update
!-------------------
ELSE IF (lev_hist==7) THEN
phi_bayes =   (1.0_dp - hk1(task7))*( (  p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)*p_success(2,task3,outcome3, priortype)&
            *p_success(2,task4,outcome4, priortype)*p_success(2,task5,outcome5, priortype)*p_success(2,task6,outcome6, priortype)&
			*p_success(2,task7,outcome7, priortype))*phi)&
           /( (  p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)*p_success(2,task3,outcome3, priortype)&
            *p_success(2,task4,outcome4, priortype)*p_success(2,task5,outcome5, priortype)*p_success(2,task6,outcome6, priortype)&
			*p_success(2,task7,outcome7, priortype))*phi&
			+p_success(1,task1,outcome1, priortype)*p_success(1,task2,outcome2, priortype)*p_success(1,task3,outcome3, priortype)&
			*p_success(1,task4,outcome4, priortype)*p_success(1,task5,outcome5, priortype)*p_success(1,task6,outcome6, priortype)&
			*p_success(1,task7,outcome7, priortype)*(1.0_dp-phi) ) + hk2(task7,outcome7,ttenure)
 

IF (phi_bayes < 0.0_dp) THEN

 phi_bayes = 1.0d-10

ELSE IF (phi_bayes > 1.0_dp) THEN

 phi_bayes = 1.0_dp - 1.0d-10

!PRINT*, "The posterior in t=8 is not interior"

END IF   


!-------------------
!h. Period-9 update
!-------------------
ELSE IF (lev_hist==8) THEN
phi_bayes = (1.0_dp - hk1(task8))*( ( p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)*p_success(2,task3,outcome3, priortype)&
            *p_success(2,task4,outcome4, priortype)*p_success(2,task5,outcome5, priortype)*p_success(2,task6,outcome6, priortype)&
			*p_success(2,task7,outcome7, priortype)*p_success(2,task8,outcome8, priortype))*phi)&
           /(  (  p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)*p_success(2,task3,outcome3, priortype)&
            *p_success(2,task4,outcome4, priortype)*p_success(2,task5,outcome5, priortype)*p_success(2,task6,outcome6, priortype)&
			*p_success(2,task7,outcome7, priortype)*p_success(2,task8,outcome8, priortype))*phi&
			+p_success(1,task1,outcome1, priortype)*p_success(1,task2,outcome2, priortype)*p_success(1,task3,outcome3, priortype)&
			*p_success(1,task4,outcome4, priortype)*p_success(1,task5,outcome5, priortype)*p_success(1,task6,outcome6, priortype)&
			*p_success(1,task7,outcome7, priortype)*p_success(1,task8,outcome8, priortype)&
			*(1.0_dp-phi) ) + hk2(task8,outcome8,ttenure)


IF (phi_bayes < 0.0_dp) THEN

 phi_bayes = 1.0d-10

ELSE IF (phi_bayes > 1.0_dp) THEN

 phi_bayes = 1.0_dp - 1.0d-10

!PRINT*, "The posterior in t=9 is not interior"

END IF    


!---------------------
!i. Period-10 update
!---------------------
ELSE IF (lev_hist==9) THEN
phi_bayes =   (1.0_dp - hk1(task9) )*(  (  p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)*p_success(2,task3,outcome3, priortype)&
            *p_success(2,task4,outcome4, priortype)*p_success(2,task5,outcome5, priortype)*p_success(2,task6,outcome6, priortype)&
			*p_success(2,task7,outcome7, priortype)*p_success(2,task8,outcome8, priortype)&
			*p_success(2,task9,outcome9, priortype))*phi)&
           /( (   p_success(2,task1,outcome1, priortype)*p_success(2,task2,outcome2, priortype)*p_success(2,task3,outcome3, priortype)&
            *p_success(2,task4,outcome4, priortype)*p_success(2,task5,outcome5, priortype)*p_success(2,task6,outcome6, priortype)&
			*p_success(2,task7,outcome7, priortype)*p_success(2,task8,outcome8, priortype)&
			*p_success(2,task9,outcome9, priortype))*phi&
			+p_success(1,task1,outcome1, priortype)*p_success(1,task2,outcome2, priortype)*p_success(1,task3,outcome3, priortype)&
			*p_success(1,task4,outcome4, priortype)*p_success(1,task5,outcome5, priortype)*p_success(1,task6,outcome6, priortype)&
			*p_success(1,task7,outcome7, priortype)*p_success(1,task8,outcome8, priortype)&
			*p_success(1,task9,outcome9, priortype)*(1.0_dp-phi)   ) + hk2(task9,outcome9,ttenure)


IF (phi_bayes < 0.0_dp) THEN

 phi_bayes = 1.0d-10

ELSE IF (phi_bayes > 1.0_dp) THEN

 phi_bayes = 1.0_dp - 1.0d-10

!PRINT*, "The posterior in t=10 is not interior"

END IF            

ELSE 
!PRINT*, "Error in computing the t-step-ahead posterior"
END IF



END FUNCTION phi_bayes



!**************************************************************
!(4) Function to compute the probability of the observed level
!**************************************************************
FUNCTION problevel(level,value,value0,value1,value2,value3,ten,lev_prev)


REAL(dp), INTENT(IN):: value,value0,value1,value2,value3
INTEGER, INTENT(IN):: level, ten, lev_prev
REAL(dp):: S0, S1, S2, S3, SUMV
REAL(dp), PARAMETER:: big = 6.00_dp, small = 1.0d-35 
REAL(dp):: problevel
REAL(dp):: taut
REAL(dp):: MAXS


!IF(ten==1) THEN
! PRINT*, 'tenure', ten
!END IF
 
!-------------------------------------
!a. Compute the component differences
!-------------------------------------
IF (level==-2 .or. level==1 .or. level==2 .or. level>=3) THEN

   S0 = DEXP((value0 - value)/tau)  !employment at firm C (the rest of the market)
   S1 = DEXP((value1 - value)/tau)
   S2 = DEXP((value2 - value)/tau)
   S3 = DEXP((value3 - value)/tau)

ELSE
 PRINT*, 'Error in computing the level probability'
 PRINT*, 'Tenure =' , ten
 PAUSE

END IF


MAXS = DMAX1(S0,S1,S2,S3)   



!---------------------------------------------------------------------------------------------------------
!b. Sum all these differences (through the if construct set a tolerance level to avoid division by zero)
!---------------------------------------------------------------------------------------------------------
  SUMV = S0 + S1 + S2 + S3
  
  IF(SUMV >= small) THEN 
   SUMV = SUMV
  ELSE
   SUMV = SUMV + small
  END IF
 



!--------------------------------------------------
!c. Compute retention and assignment probabilities
!---------------------------------------------------
    IF (level==-2) THEN
 
      problevel =   (1.0_dp - seprate(lev_prev,ten-1))*S0/SUMV+ seprate(lev_prev,ten-1) 

      
    ELSE IF (level==1) THEN


      problevel =   (1.0_dp - seprate(lev_prev,ten-1))*S1/SUMV
	   


    ELSE IF (level==2) THEN

    

      problevel =    (1.0_dp - seprate(lev_prev,ten-1))*S2/SUMV

     

    ELSE IF (level>=3) THEN

    
     problevel =   (1.0_dp - seprate(lev_prev,ten-1))*S3/SUMV
	

    END IF 


   ! correction for low probability events
   IF(problevel < small) THEN
     problevel = small
   END IF


END FUNCTION problevel





!******************************************************
!(5) Function for the separation rate at each level
!******************************************************
FUNCTION seprate(level, tensep)

INTEGER, INTENT(IN):: level, tensep
REAL(dp):: seprate


IF(level==-2 )THEN
 seprate = 0.0_dp
 PRINT*, 'Error with separation shocks (past level -2)'
 PAUSE

!-----------
!a. Level 1 
!-----------
ELSE IF(level==-4)THEN
 seprate = 0.0_dp
 !PRINT*, 'error with separation'
 
ElSE IF(level==1 .and. tensep==0)THEN
 seprate = 0.0_dp

!xi_21
ElSE IF(level==1 .and. tensep==1)THEN
 seprate = xi_1

ElSE IF(level==1 .and. tensep==2)THEN
 seprate = xi_1

ElSE IF(level==1 .and. tensep==3)THEN
 seprate = xi_11 + plus

ElSE IF(level==1 .and. tensep==4)THEN
 seprate = xi_11

ElSE IF(level==1 .and. tensep==5)THEN
 seprate = xi_21

ElSE IF(level==1 .and. tensep==6)THEN
 seprate = xi_21

ElSE IF(level==1 .and. tensep==7)THEN
 seprate = xi_21

ElSE IF(level==1 .and. tensep>=8)THEN
 seprate = xi_21


!-----------
!b. Level 2 
!-----------
ElSE IF(level==2 .and. tensep==0)THEN
 seprate = 0.0_dp


ELSE IF(level==2  .and. tensep == 1)THEN
 seprate = xi_2

ELSE IF(level==2  .and. tensep == 2)THEN
 seprate = xi_2

ElSE IF(level==2 .and. tensep==3)THEN
 seprate = xi_2 + plus

ElSE IF(level==2 .and. tensep==4)THEN
 seprate = xi_14

ELSE IF(level==2  .and. tensep == 5)THEN
 seprate = xi_34

ELSE IF(level==2  .and. tensep == 6)THEN
 seprate = cpen


ELSE IF(level==2  .and. tensep == 7)THEN
 seprate = xi_25

ELSE IF(level==2  .and. tensep >= 8)THEN
 seprate = xi_25


!------------
!c. Level 3 
!------------
ElSE IF(level>=3 .and. tensep==0)THEN
 seprate = 0.0_dp


ELSE IF(level>=3  .and. tensep == 1)THEN
 seprate = xi_3

ELSE IF(level>=3  .and. tensep == 2)THEN
 seprate = xi_3

ElSE IF(level>=3 .and. tensep==3)THEN
 seprate = xi_3

ElSE IF(level>=3 .and. tensep==4)THEN
 
 seprate = xi_14 

ELSE IF(level>=3  .and. tensep == 5)THEN
 
 seprate = xi_34 

ELSE IF(level>=3  .and. tensep == 6)THEN
 seprate = cpen 

ELSE IF(level>=3  .and. tensep == 7)THEN
 seprate = xi_25 


ELSE IF(level>=3  .and. tensep >= 8)THEN
 seprate = 0.0_dp 

ELSE 
 PRINT*, "Error in computing the separation rate (not a level observed)"
END IF


END FUNCTION seprate




END MODULE mod_update